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A method of modelling the three-dimensional microstructure of random isotropic two-phase ma- 
terials is proposed. The information required to implement the technique can be obtained from 
two-dimensional images of the microstructure. The reconstructed models share two-point corre- 
lation and chord-distribution functions with the original composite. The method is designed to 
produce models for computationally and theoretically predicting the effective macroscopic prop- 
erties of random materials (such as electrical and thermal conductivity, permeability and elastic 
moduli). To test the method we reconstruct the morphology and predict the conductivity of the 
well known overlapping sphere model. The results are in very good agreement with data for the 
original model. 

47.55. Mh, 44.30.-Hv, 81.05. Rm, 61.43. Bn 



Predicting the macroscopic properties of composite 
or porous materials with random niicrostructures is an 
important problem in a range of fields |h,p[. There now 
exist large-scale computational methods for calculating 
the properties of composites given a digital representa- 
tion of their microstructure (eg. permeability ^^ , con- 
ductivity [p|-pi and elastic moduli 0). A critical prob- 
lem is actually obtaining an accurate three-dimensional 
description of this microstructure P,nJ8| . For particular 
materials it may be possible to simulate microstructure 
formation from first principles. Generally this relies on 
a detailed knowledge of the physics and chemistry of the 
system; the accurate modelling of each material requiring 
a significant amount of research. Where such information 
is unavailable an alternative is to directly |g|-|l5|| or statis- 
tically [pl0,@,Eq-El| reconstruct the microstructure from 
experimental images. 

Several techniques of direct reconstruction have been 
implemented. A composite can be repeatedly sectioned, 
imaged and the results combined to reproduce a three- 
dimensional digital image of the microstructure [pHll|. 
For porous materials, time consuming sectioning can 
been avoided by using laser microscopy |l2| which can 
image pores to depths of around 150^m. Recent micro- 
tomography studies have also directly imaged the three- 
dimensional microstructure of porous sandstones [ [l3|Jl4| 
and magnetic gels |15| . The complexity and restrictions 
of these methods provide the impetus to study alterna- 
tive reconstruction methods. 

Based on the work of Joshi ||l^, Quiblier [|l^ intro- 
duced a method of generating a three-dimensional statis- 
tical reconstruction of a random composite. The method 
is based on matching statistical properties of a three- 
dimensional model to those of a real microstructure. A 
key advantage of this approach is that the required in- 
formation can be obtained from a two dimensional im- 
age of the sample. Recently the method has been ap- 
plied to the reconstruction of sandstone P,p|,p^|jl9|] and 



a material composed of overlapping spheres Q. Com 
putations of the permeability and conductivity gl 18 
of the reconstructed images underestimate experimental 
data by around a factor of three. This can be partially 
attributed to the fact that percolation threshold of the 
reconstructed models is around 10% while the experimen- 
tal systems had thresholds of less than 3% Q. Recent 
work in microstructure modelling has led to a general 
scheme [p|,p2|-p7t (§ |) which includes the model employed 
by Quiblier. Importantly, other models in the scheme can 
mimic the low percolation thresholds observed in sand- 
stones (and many other materials ||2^). It is therefore 
timely to reconsider statistical methods of reconstruct- 
ing composite microstructure. 

Prior methods of statistical reconstruction produce 
three-dimensional models which share first (volume frac- 
tion) and second (two-point correlation function) order 
statistics with the original sample. However the com- 
plete statistical description of a random disordered ma- 
terial requires higher order information [g|,^ (eg. the 
three and four point correlation functions). Information 
which in turn is a crucial ingredient of rigourous theories 
of macroscopic properties ||l|,^,^ , and therefore impor- 
tant to the success of the model. In this paper we show 
that reconstructions based on matching first and second 
order statistics do not necessarily provide good models of 
the original composite (§ O). An alternative method of 
reconstruction is proposed and tested (§ III). The pro- 



cedure is employed to reconstruct a composite generated 
from identical overlapping spheres (lOS) and successfully 
predicts the electrical conductivity of the model (§ [rv|) . 



I. MODEL COMPOSITE MATERIALS 



To study the statistical properties of composites it is 
conventional to introduce a phase function 4>{r) which 
equals unity or zero as r is in phase one or two. The 
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volume fraction of phase one is p — (0), while the stan- 
dard two-point correlation function is defined as p"' (r) = 
{(j){ri)(p{r2)) with r=\r2 — f-ij (assuming the material is 
statistically homogeneous and isotropic) . p^'^' (r) repre- 
sents the probability that two points a distance r apart 
will lie in phase one. From the definition p"-* (0) = p and 
linir^oo P^'^' (r) — *■ p^. The surface area per unit volume 
is s = — 4(ip(^^/dr|,.=o @1- Higher order functions can 
be analogously defined, these playing a central role in 
rigourous theories of composite properties P8| . In prac- 
tice the correlation functions of real composites beyond 
second order are difficult to measure and there are signifi- 
cant advantages in developing models for which the func- 
tions are exactly known. The primary models in this class 
are the identical overlapping sphere model (lOS) |^l|], its 
generalisation to overlapping annuli p^J and models de- 
rived from Gaussian random fields (GRFs) ||5|, p2| , p2| , p3[ 
which are central to reconstruction procedures. 

We utilise two methods of generating isotropic GRFs. 
Each has specific advantages which we discuss. The first 
method develops the random field in a cube of side- length 
T using a Fourier summation; 



terms of the first definition P{k) = Ank? p{k). In this 
case the fields are not periodic, but N can be chosen 
arbitrarily largely over a specified k range. This is espe- 
cially useful for resolving p{k) (so that Eq. (||) holds) in 
cases where it is strongly spiked (eg. P{k) = 5{k)) [ p3[ . 

Model composite materials can be defined from a GRF 
y{r) by taking the region in space where a < y{r) < j3 
as phase one and the the remaining regions {y{r) < a 
and y(r) > (3) as phase two. This is the "two-level cut" 
random field of Berk |Q. In the case a = — oo the 
more common "one- level cut" field is recovered P,0,K2i . 
The phase function of this model is (/)(y(r)) — H{y{r) — 
a) — H{y{r) — (3) where H is the Heaviside step func- 
tion. The joint probability distribution of the corre- 
lated random variables y = [j/(ri), ?;(r2), . . . , y(r„)]-^ is, 
P„(y) = ((27r)"|G'|)-2exp(-iy^G~it/) where the ele- 
ments of G are g^ = gi^ij) — {y{'^i)y{'^j))- Therefore the 
n-point correlation function is 
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where ki. 



^imn = ^{li + mj + nk). The statistics 
of the field are determined by the random variables 
ci,m,n=aimn + ihrnn {.dimn and 5/m„ real). We require 
that y is real (Q,m,„=c-;, _„,_„) and that {y) = 
(co, 0,0=0). To ensure isotropy we also take cimn=0 
for fcim„ = |A;im„|>27riV/T. To generate a Gaussian 
field the coefficients «;„„ are taken as random inde- 
pendent variables (subject to the conditions on cimn) 
with Gaussian distributions such that (aimn) — and 

(aLn)=5P(^imn) (t") (similarly for kmn)- The func- 
tion p{k) is a spectral density. It can be shown that a 
random field defined in this manner has field-field corre- 
lation function 



g{r) = (y{ri)yir2)) 
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By convention g(0) = 1 which sets a constant of pro- 
portionality on p{k). The definition (m) can be efficiently 
evaluated by an FFT routine ^ and is T-periodic in each 
direction. This is valuable for approximating an infinite 
medium in calculations of macroscopic properties. 

Alternatively a random field can be generated using 
the "random-wave" form |3lJ32l 



The volume fraction of phase one is p = p^^' — h = {pp — 
Pa) where pa — (27r)~2 J"^ e~* /"^dt and p^^^(r) = h{r) 
with |32l|3|] 
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The auxiliary variables h and h{r) are needed below. The 
three-point correlation functions |28] have also been eval- 
uated IJH. 

We now show how new models can be developed. Sup- 
pose (t>i{r) and 4)2{i') are the phase functions of two 
statistically independent composites with volume frac- 
tions pi and p2 and correlation functions p\ and pj ■ 
New model composites can be formed from the inter- 
section and union sets of each structure. The inter- 
section set (/)(r) = 01 (r) x (t>2{r) has volume fraction 
p= {4'i(r)4'2{r)) = {4>i{r)){(i)2{r)) =PiP2 and correlation 
function 

p(2)(r) ^ (Mri)Mri)Mr2)Mr2)) (6) 

= (0i(ri)(/)i(r2))(02(ri)02(r-2)) = pf\r)p^2\r). 
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where (pi is a uniform deviate on [0, 27r) and ki is uni- 
formly distributed on a unit sphere. The magnitude of 
the wave vectors ki are distributed on [0, oo) with a prob- 
ability (spectral) density P{k) {J^ P{k)dk = 1). In 



In a similar way a composite can be modelled as the 
union of two independent models. In this case the phase 
function is (f>{r) — (f>i{r) + 02 (r) — (j)i{r)(j)2{r) so that 
p^Pi+P2- P1P2 and 



^(^\r) ^ pf\r){l - 2p2) +pf (r)(l - 2p,) 



+2pip2+pl'^(r-)pf (r). 
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Therefore if the statistical properties of the original mor- 
phologies are known (eg. level-cut GRF's or the over- 
lapping sphere model) the properties of their union and 
intersection sets are also known [g7|. Note that these 
results apply to arbitrary independent phase functions, 
are simply extended to three or more independent sets, 
as well as to the calculation of higher order correlation 
functions. These simple results greatly extend the classes 
of morphology which can be reproduced by the models. 

To simplify matters we now restrict attention to a few 
primary models of microstructure. Consider first struc- 
tures derived using the normal two-level cut GRF scheme 
(model "N" ) . These have the basic statistical properties 
p—h (recall h=pi3 ~ Pa), p^'^\r) — hir) and s — — 4/i'(0). 
We also take pa = c x (1 — p)/2 and pp — Pa + P 
(c € [0,1]) to specify the level-cut parameters; for ex- 
ample c = corresponds to a one-cut field {pa,0 = 0,p 
or a = — oo) and c = 1 to a symmetric two-cut field 
{Pa,i3 — i~f'5 + f or a — — /3). Second, we take a 
class of models based on the intersection set (model "I" ) 
of two statistically identical level-cut GRF's. For this 
model p = h?, p'^'^\r) = h'^{r) and s — — 8^/i'(0) with 
Pa = c(l — s/p)!^ and Pp = Pa + y^. Finally, we intro- 
duce a model based on the union set (model "U" ) of two 
two level-cut fields. In this case p 

2h^ + 2 h{r){l - 2h) + K^jr) and s -- 

Pa = c^l -p/2 and p^g = Pa + 1 - VI ^ P- 

To generate examples of the models defined above we 
employ the field-field correlation function p7^,p3,p6| 



-8VT^/i'(0) with 
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characterised by a correlation length ^, domain scale d 
and a cut-off scale rr ■ This has Fourier transform 



p(fc) = 



\i~r,)-^ed^ 



[d^ + e{kd - 27r)2][d2 -K i;^{kd + 27r)2 



(9) 
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[d2 + rl{kd - 27r)2][d2 + rl{kd + 27r)2] ' 

Note that g{r) is symmetric in Vc and ^ and remains well 
defined in the limits r^ — > f and r^ or ^ — > c». In the lat- 
ter cases p{k) -^ S{k — 2TT/d)/Ank^ |33|. For the purposes 
of calculating the surface area. 
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In the case re or ^ ^ a fractal surface results p5| , |33[| 
Cross-sections of six of the model microstructures ob- 
tained with Tc—l, ^—2 and d—2iim are illustrated in 
Fig. y. p^2)(|.) ig measured from three-dimensional re- 
alisations (using 128'^ pixels) of the models and plotted 
against its theoretical value in Fig. 0. The agreement is 
very good. In the following section we also consider each 
of the models at an intermediate value of c = 2- The 
extra three models, along with the six shown in Fig. |l] 




(a) N,c=0, p=0.4 (b) I,c=0,p=0.16 (c) U,f,=0, p=0.64 








((i) N,c=l,p=0.4 (e) I,c=l,p=0.16 (f) U,r,= l, p=0.64 

FIG. 1. Six different microstructures generated by the 
level-cut scheme. In the top row we show a one-cut field 
and its intersection and union with a statistically identical 
structure. In the bottom row we show analogous structures 
derived from a two-cut field. The images have a side length 
of lO/xm. 




1.5 2.0 2.5 3.0 
r(|im) 

FIG. 2. The theoretical (lines) and measured (symbols) 
correlation functions of the six models shown in Fig. hi The 
squares correspond to the models constructed from one-cut 
fields (Figs, tk-c), and the triangles to the the two-cut fields 
(Figs.jild-f). 

give nine primary classes of microstructure with which to 
compare real composites. These broadly cover the types 
of morphology obtainable by combining two composites 
generated by the level-cut GRF scheme. 

II. STATISTICAL RECONSTRUCTION 

The two most common experimentally measured mor- 
phological quantities of composites are the volume frac- 

(2) 

tion Pexp and the two-point correlation function pexp{r) 
(eg. Refs. |,|l|,^,|^,|8|). Consider how this information 
might be used to reconstruct the composite using the 
simple one-cut GRF model (model N, c = or a = —oo). 
The level-cut parameter P can be obtained by solving 
Pexp = (27r)~2 J_ e~* '"^dt and the field-field function 
obtained by numerical inversion of 



p(2) (r) =p^ + 
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From g{r) we can obtain p{k) by inverting Eq. (0) and 
using either Eq. (|l|) or (||) to obtain y(r) and hence the 
model phase function 0(r). The reconstruction shares 
first and second order statistical properties with the im- 
age and would therefore be expected to yield a reason- 
able model of the original composite. This is similar to 
the procedure of Quiblier p^ employed in previous stud- 
ies p,[I|Jq,p!q-pi|| although the formulation of the model is 
different. There are several operational problems with 
this reconstruction procedure. First, the numerical in- 
version of Eq. (O) may not be robust or well defined. 

— (2) 

Furthermore experimental error in pexpir) is carried over 
to g{r). Second, the inversion of Eq. (11) may yield a 
spectral density p{k) which is not strictly positive. We 
now generalise the method to incorporate the models N, I 
and U of § y and show how these problems can be avoided. 
First select one of the three models (N, I or U) and a 
value of c=0, ^ or 1 (giving a total of nine combinations) 
so that a and (3 are fixed by Pexp- It remains to find g{r). 
Instead of inverting an analog of Eq. ( |ll| ) we assume this 
function is of the general form given by Eq. (g) (this 
guarantees that p{k) is positive). The three length scale 
parameters are obtained by a best fit procedure which 
minimises the normalised least-squares error; 
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Here p\^^{ri) = p^ > [g{ri] r^ ^, d)] is the correlation func- 
tion appropriate for model N, I or U. Once r^, S, and d 
have been obtained the reconstruction 4>{r) can be gener- 
ated. If the one-cut model (N, c = 0) is chosen we assume 
that the results will not differ significantly from those 
obtained using Quiblier's method. 

To illustrate the procedure we reconstruct a material 
with known statistical properties. For this purpose we 
choose a normal two-cut GRF model with p^.p — 0.4, 0.6 
(ie. model N, c = 1) obtained from the field-field func- 
tion § 
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with Iq = 2.0^m. The "experimental" data for the recon- 
struction pexp{ri) are evaluated using Eq. (||) at 80 points 
distributed uniformly on the interval [0,4]/im (shown as 
symbols in Fig. H). The minimisation algorithm is used 
to find Tc, ^ and d for four different models. Numerical 
results are reported in Table I and the best-fit functions 
pJ^ are plotted in Fig. ||. Each of the models is able to 
provide an excellent fit of the data. As expected, model 
N (c = 1) provides the least value of Ep^^\ However the 



TABLE I. Theparameters obtained in the reconstruction 
procedure (Eq. (121)) of a test composite. The surface area 
of the original model is O.S7pui~^ . Here, and in subsequent 
tables, n{m) denotes nxlO"'". 
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Experimental data 
Model N (c=0) 
Model N (0=1) 
Model I (c=l) 
Model U (c=l) - 
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FIG. 3. The correlation functions pL^(r) (lines) of four re- 
constructed models obtained by fitting "experimental" data 
(symbols) . 

relative improvement over the other three models is not 
large, and probably of little significance in the presence 
of experimental error. Cross-sections of the original com- 
posite and the reconstructions are shown in Fig. H(a)-(e). 
The extremely different morphologies exhibited by the re- 
constructions provide a graphical illustration of the non- 
uniqueness of p'^2) (j-) ^ Therefore for prediction of macro- 
scopic properties (which will differ dramatically for ma- 
terials shown in Fig. 0) it is necessary to find a more dis- 
criminating method of distinguishing composites. From 
the cross-sectional images the best candidates appear to 
be models N (c= 1) and U (c= 1) shown in Figs, ^(c) and 
(e). Obviously it is preferable to establish some quanti- 
tative test to choose the best representation. 

A second useful illustration of the method is pro- 
vided by reconstructing a material with a strongly os- 
cillating correlation function. For this case we take as 
a test-composite a one-cut model with p = 0.2 and 
Pa,i3 — 0.0^0.2 (ie. model N, c = 0) based on the field-field 
function 



^(r) = 3r ^{k\~kl) ^(sin fcir - sinfcor) 



— 3r 2(fc3_fc3^ ^(fci cos/cir — fcocosfco'") 
p{k) = 3[47r(fc3 _ kl)]-^[H{k - fco) - H{k - fci)] 



(14) 
(15) 

with fco = 3.0 and fci = 4.5(/im)"^. The oscillatory be- 
haviour of the correlation function (see Fig. g) can only 
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(a) Exp. 



(b) N, f-0 (c.) N, r-1 (d) I, r-1 (e) U, f-1 




{f)Bxp. {g)I,c = 0{h)I,c = 0(i)U,r. = Oa)N,r,= | 

FIG. 4. Realisations of "experimental" and reconstructed 
composites. Top row: A material with a monotonicly decay- 
ing correlation function (a) compared with four reconstruc- 
tions (b)-(e). The two point correlation functions of each 
composite are practically identical (see Fig. pi). Bottom row: 
A model composite exhibiting an oscillatory correlation func- 
tion (f) and four reconstructions (g)-(j). In each case the 
region shown is lOxlO/im. 
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FIG. 5. Correlation functions of two reconstructions (lines) 
of a material exhibiting an oscillatory p'^-'(r) (symbols). A 
"mild" two-cut model (dashed line) is unable to accurately 
reproduce the strong oscillations. 

be reproduced by three of the nine basic microstructures; 
models N,I and U with c — (ie. those formed from 
one-cut fields). For these models Ep'-^^ < 0.005 whereas 
Ep'^^^ >0.02 for those based on two-cut structures (c> ^ 
so < Pa < p/i). To illustrate this we show the best 
fit of a normal two-cut model with pa^p = 0.05, 0.25 (N, 
c= |). As can be seen in Fig. || this "mild" two-cut model 
(shown as a dashed line) cannot reproduce the behaviour 
of the experimental data (see Table O) . Realisations of 
the original material and reconstructions are shown in 
Fig. H(f)-(j). Each appears to provide a reasonable rep- 
resentation. 

In contrast to the case of a monotonicly decaying 
P {r) (which was reproduced by four distinct models) 
strong oscillations appear to be a signature of morpholo- 
gies generated by the single level-cut model. Unless there 
exists some reason to employ models U and I in such a 
case it is likely that the standard one-cut GRF (ie. the 
model employed in prior studies) will be appropriate. 



TABLE II. Reconstruction of a normal one-cut model 
with an oscillatory correlation function. Models formed from 
two-cut fields (ie. pa > 0) were unable to reproduce the os- 
cillations of p"'(r) (see eg. row 4). The surface area of the 
original model is l.OO/im"^. 
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There is also a physical basis for this argument when 
spinodal decomposition plays a role in the microstruc- 
tural formation. In this case Cahn |^ has shown that 
the evolution of the phase interface is described by the 
level-set of a sum of random waves similar to (||) . 

Finally we comment on the morphological origin of the 
oscillations, and why they cannot be well reproduced by 
two-cut models. In Fig. we show p^ \r) and an image 
of model N, c = with rc=2, ^=4 and d—lfixa.. The 
material has strong oscillatory correlations, these repre- 
senting the "regular" alternating domains which appear 
in the image. Compare this with data shown for the two- 
cut model (N, c= i) obtained from the same GRF: the 
alternating structure is still present but the oscillations 
are practically extinguished. This is due to the sharper 
decay (or equivalently the doubled specific surface) as- 
sociated with the thinner two-cut structures [M. For 
comparison we also show a structure with no repeat scale 
(model N, c = with rc=\, ^=\ and d=100/xm). 
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FIG. 6. Three different types of microstructure. A one-cut 
model with a well defined domain (or repeat) scale (left), a 
two-cut model obtained from the same GRF (centre) and a 
one-cut field with no domain scale (right). The oscillations 
of p' ' are very weak for the central model even though the 
domain scale is obvious to the eye. 
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III. COMPARISON OF HIGHER-ORDER 
STATISTICAL PROPERTIES 



We have shown that reconstructions exhibiting quite 
different morphological properties can share the same 
two-point correlation function. Here we propose and test 
three methods with the aim of finding a way of select- 
ing the best reconstruction. Following Yao et al H we 
can compare the three-point correlation function of the 
model and experimental materials. To do so we define a 
normalised least-square measure of the error as 

Nr N, Ng 

i=l j=l k=l 

Nr N^ Ng 

- E E Ei^-U^- ^^■' ^fc) - p-p]'- (16) 

i=l i=l fc=l 

The three-point function p^^'{r, s,9) gives the probabil- 
ity that three points distances r, s and t — {r'^ + s'^ — 
2rscos6)^ apart all lie in phase one. For our examples 
we take Nr.s,e = 8 with a uniform distribution of r and 
s on [0, 2]/xm and 9 on [0, n] 

A second method of characterising morphology is to 
calculate microstructure parameters which appear in the- 
oretical bounds on transport and elastic properties ||l|,^ . 
We therefore expect the parameters to contain critical in- 
formation about the aspects of microstructure relevant to 
macroscopic properties. These are 



but it can be simply evaluated from realisations of the 
experimental and reconstructed materials. To quantify 
the difference between the chord-distributions we again 
employ a least-squares error; 
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where f{r,s,9) — p^'^^ (r, s, 9) —p^^^ {r)p^^^ {s)/p, q — 1 — P, 
u — cos 9 and Pn{u) denotes the Legendre polynomial of 
order n. The parameter C, occurs in bounds on the con- 
ductivity and the bulk modulus, while rj occurs in bounds 
on the shear moduli. As pJ^ and peJp are available for 
our test models the parameters can be calculated P,p2[. 
Techniques have also been suggested for directly evalu- 
ating the parameters from experimental images pO| , pd| . 
We anticipate that the closer C^fn is to C,f^xp the better 
the reconstructed model. Note that C, and 77 contain only 
third order statistical information and higher order in- 
formation is potentially important for our purposes. 

A third simple measure of microstructure is the chord- 
distribution function of each phase 1 40 4^:^§1 . For phase 
one this is obtained by placing lines through the com- 
posite and counting the number of chords n{r) of a given 
length r which lie in phase one. The chord-distribution is 
defined as p'''^\r) = n{r)/ J^ n{r)dr so that p'^^\r)dr is 
the probability that a randomly selected chord will have 

length between r and r + dr. p^'^\r) is defined in an anal- "Exp " data 0^2 OM (Tst" 
ogous manner. At present it is not possible to analyti- 
cally evaluate this function for the level-cut GRF media, 



with j = 1,2. Note that p'^^\r) contains information 
about the degree of connectedness in phase j and thus 
is likely to incorporate important information regarding 
macroscopic properties Q . 

We also compute the conductivity of samples (size 128'^ 
pixels) using a finite-difference scheme g] . We choose the 
conductivity of phase one as txi = 1 (arbitrary units) and 
phase two insulating ((T2 = 0). At this contrast the ef- 
fective conductivity a is very sensitive to microstructure. 
The results therefore allow us to gauge the ability of a 
reconstruction to predict macroscopic properties. This 
contrast also occurs commonly in a range of systems (eg. 
electrical conductivity of brine saturated porous rocks or 
thermal conductivity of aerogels and foams). 

We have calculated the morphological quantities de- 
fined above for the first four reconstructio ns ( reported in 
Table |). The results are shown in Table III. First note 
that i?p(^) is greater than Ep'^^^ by a factor of 2-5 ^] 
in each case and is probably of little use in an actual 
reconstruction. The values of the microstructure param- 
eters C and 77 are conclusive; as we expect they indicate 
that model N (c = 1) is best. The chord-distributions of 
the experimental and reconstructed material are shown 
in Fig. (phase one) and Fig. ^ (phase two). From Ta- 
ble |l| we see that the chord-distribution provides a very 
strong signature of microstructure. The results indicate 
that either model N (c= 1) or model U (c= 1) is the best 
reconstruction. The fact that the conductivity of each 
model is so close to the experimental data provides some 
evidence that matching the chord distributions is more 
important than matching C, and 77. The same compari- 
son is shown for the reconstructions of the test composite 
which exhibits an oscillatory p^^^ (r) in Table IV . Model 
N (c = 0) provides the best reconstruction based on both 
the chord-distribution and the microstructure parame- 
ters. This leads to a good prediction of the conductivity. 

TABLE III. A comparison of the statistical and transport 
properties of the four reconstructed models (Table M with 
those of the "experimental" composite. The measured sur- 
face area of the digital reconstructions is also shown. 



CI 



Ep 



(3) 



C/i: 



■nfit 



Ep 



(1) 



Ep^ > CTree /o"l 



N 





5(-3) 


0.32 


0.29 


1.06 


0.25 


0.62 


0.032 


N 


1 


9(-5) 


0.74 


0.54 


0.75 


0.04 


0.11 


0.114 


I 


1 


2(-3) 


0.47 


0.37 


0.98 


0.20 


0.48 


0.069 


U 


1 


6(-3) 


0.87 


0.70 


1.02 


0.02 


0.15 


0.120 



0.110 



Exp. data 
Model N (c=0) 
Model N(c=l) 
ModelI(c=l) 
Model U (c=l) 




1.5 2.0 2.5 3.0 
r(|im) 

FIG. 7. The chord-distribution (for phase 1) of an "experi- 
mental" composite (Fig. ^) compared with data for the four 
reconstructions shown in Figs. 0(b)-(e). Both models N and 
\J {c — 1) appear to mimic the "experimental" data. The 
lines in the graph are guides to the eye only. 



Tables [11 and IV, p^ ' (like p^ ') does not appear to pro- 
vide a strong signature of microstructure [Q . It is there- 
fore not possible to conclude that a good reproduction of 
p(3) (qj. j)(4)) implies a successful reconstruction as was 
done in Rcf. m . In contrast both the chord-distributions 
and the microstructure parameters appear to provide a 
strong signature of composite morphology, and hence a 
method of selecting a useful reconstruction of the original 
material. 

TABLE IV. A comparison of the statistical and transport 
properties of the three reconstructed models (Table |l|) which 
are able to reproduce the oscillatory correlation function of 
a test composite. 
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i5p(3) 


C/.t 


Vfit 


Srec 


i5p(i) 


£p(2) 


CTrec/cTl 


N 
I 
U 


9(-5) 
5(-3) 
5(-3) 


0.24 
0.33 
0.20 


0.20 
0.25 
0.17 


1.00 
1.16 
1.10 


0.001 
0.137 
0.008 


0.003 
0.036 
0.127 


0.025 
0.032 
0.009 


"Exp." 


data 
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0.023 
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IV. RECONSTRUCTION OF THE lOS MODEL 
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FIG. 8. The chord-distribution (for phase 2) of an "exper- 
imental' composite compared with data for four reconstruc- 
tions (see caption of Fig. m). 

In § P we showed that it was possible to generate a 
number of morphologically distinct reconstructions which 
share first and second order statistical properties with an 
experimental composite. Here we have suggested three 
methods of choosing the best reconstruction. As Ep^^^^ 
is relatively small for all seven reconstructions shown in 



Realisations of the identical overla ppi ng sphere (lOS) 
model |3^] (or Poisson grain model |46|) are generated 
by randomly placing spheres into a solid or void. In 
the latter case the morphology is thought to provide 
a reasonable model of the pore-space in granular rocks 
(so transport occurs in the irregular void region). As 
the model has a different structure to the level-cut GRF 
model it provides a useful test of reconstruction proce- 
dures ^ . The correlation function of the material [|Tj is 
p(2) (^) — pv(r) fQj. J. ^ 2ro and p^^^ (r) = p'^ for r > 2ro 
where 



^^^) = ^n(^)^^fe 



(20) 



For this model it is also possible to calculate the pore 
chord distribution as p^^\r) ^ -3/4ro x lnj5p3r/4ro Q^ 
We first consider the lOS model at volume fraction 
Pexp = 0.2. The system is 80% filled with spheres 
of radius ro=l/im. Nine reconstructions are generated 
(by minimising Ep^"^'), and their higher order statisti- 
cal properties are compared with those of the lOS model 
in Table Based on Ep'-^^ (and Ep'-^^) we note that 
models U (c = ^,1) perform poorly while the standard 
one-cut model is very good. The microstructure param- 
eters C and 1] indicate that the best reconstruction is 
model I (c = 1) followed by model I (c = i). However 
both models fail to reproduce the solid chord distribu- 
tion {Ep'''^^ > 0.6) which is better mimicked by models I 
(c = 0) and N (c = 0) . The ambiguity of the results indi- 
cate that none of models considered may be appropriate. 



TABLE V. A comparison of the statistical properties of 
11 reconstructions with those of the lOS model at porosity 
20%. Most of the models are able to reproduce the low order 
statistical properties of the lOS model. 
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Ep 
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£^P'^' s/.t C/»t r)fu Ep'-^^ Ep- 
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N 





l(-4) 


9(-4) 


0.94 


0.31 


0.28 


0.066 


0.26 


N 


1 

2 


3(-3) 


5(-3) 


0.79 


0.74 


0.54 


0.35 


0.15 


N 


1 


2(-3) 


8(-3) 


0.79 


0.84 


0.63 


0.59 


0.31 


I 





2(-4) 


7(-4) 


0.98 


0.35 


0.30 


0.024 


0.24 


I 


1 

2 


6(-4) 


l(-3) 


1.07 


0.50 


0.38 


0.042 


0.65 


I 


1 


4(-4) 


l(-3) 


1.05 


0.52 


0.40 


0.030 


0.63 


U 





2(-4) 


l(-3) 


0.92 


0.28 


0.26 


0.077 


0.30 


U 


1 

2 


l(-2) 


2(-2) 


0.91 


0.79 


0.62 


0.49 


0.11 


u 


1 


l(-2) 


2(-2) 


0.91 


0.87 


0.70 


0.40 


0.15 


Is 




7(-4) 


6(-4) 


1.00 


0.40 


0.33 


0.003 


0.23 


Iio 




l(-3) 


5(-4) 


1.00 


0.43 


0.35 


0.003 


0.13 



"Exp." data (lOS) 



0.96 0.52 0.42 



The lOS model can be thought of as the intersection 
set of infinitely many composites comprised of a single 
sphere of phase 2 (so 4>{r) = within the sphere). This 
suggests that the morphology may be better modelled 
with the level-cut scheme by increasing the number of 
primary composites beyond two. To this end we gener- 
alise model I to the case of n independent one-cut fields 
so that p'-^\r) = h"{r) with pa = 0, pp — p(^/") and 
s := ~4np^-^/"h'{0). This is termed model "I„". The 
statistical properties of the reconstructions for the cases 
n=5 and n=10 are shown in rows 10 and 11 of Table M. 
The models reproduce the "experimental" pore chord dis- 
tribution very well, and offer a progressively better repre- 
sentation of the solid chord distribution and microstruc- 
ture parameters. The chord-distributions of model I5 are 
shown in Fig. |^ along side those of the standard one-cut 
model and the lOS model. The good agreement between 
the measured and theoretical value of p^^-* (r) for the lOS 
model demonstrates the accuracy with which this func- 
tion can be evaluated for a sample of 128^ pixels. 

To determine which morphological measure ((^ and rj 
or Ep^^'> and Ep'^'^^) should be used to select the best re- 
construction we examine the model morphology and con- 
ductivity. Three-dimensional images of models N (c = 0), 
I (c=l) and /lo are shown alongside the lOS model in 
Fig. |l^. The pore space of the single-cut GRF (Fig. [ic|a) 
is more disconnected than that of the lOS model, while 
the pores are too large and uniform in the intersection 
model (Fig. |l^). Model Iio (Fig. |l^c) appears better 
able to reproduce the interconnected structures charac- 
teristic of overlapping spheres. The results for the con- 
ductivity are, ct = 0.038 for model N (c = 0), a = 0.080 for 
model I (c=l), (7 = 0.052 for model Iio and cr = 0.063 for 
lOS. The fact that model Iio better mimics lOS morphol- 
ogy and conductivity than model I (c= 1) provides evi- 
dence that minimising Ep^^' should be given more weight 
than matching experimental values of C and 77. 




2 3 

r(|im) 

FIG. 9. The chord distribution of the lOS model (open 
symbols), model Iio (solid symbols) and the standard one-cut 
model (broken line, symbols omitted for clarity). The heavy 
line is the theoretical curve for the lOS model and the lighter 
lines are guides to the eye only. 
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N^i 




(c) Recon. I 



(d) lOS Model 



FIG. 10. Reconstructions of the overlapping sphere (lOS) 
model at porosity p = 0.2. To aid visualisation the pores are 
shown as solid, and solid as void. The images shown here and 
the chord distributions (Fig. 0) indicate model Iio provides 
the best reconstruction of the lOS model. 



TABLE VI. The results of the reconstruction procedure 
for the lOS model. The specific surface of the lOS model 
is s = 0.71, 0.96, 1.08, 1.10/im~ as p increases. Generally 
model Iio provides a better match of the chord distributions 
than model I5. In each case Pi3~p^ for model I„. 



P 



CI 



^ 



d 



Ep 



(2) 



^/it 



Ep''^ Ep' 



(2) 



0.1 



0.2 



0.3 



0.4 



I5 


0.8770 


0.8769 


3.8336 


3(-4) 


0.69 


0.011 


0.33 


Iio 


1.2472 


1.2470 


3.8608 


5(-3) 


0.70 


0.011 


0.31 


I5 


0.9942 


0.9947 


3.9055 


8(-4) 


1.00 


0.003 


0.23 


Iio 


1.4173 


1.4174 


3.9777 


l(-3) 


1.00 


0.003 


0.13 


I5 


1.0974 


1.0973 


3.9756 


l(-3) 


1.14 


0.003 


0.23 


Iio 


1.6047 


1.6053 


4.0375 


l(-3) 


1.13 


0.003 


0.19 


Is 


1.2148 


1.2151 


4.0250 


l(-3) 


1.17 


0.006 


0.16 


Iio 


1.8146 


1.8158 


4.1244 


l(-3) 


1.15 


0.004 


0.18 



We adopt this strategy to reconstruct the lOS model 
at p — Q.l, 0.3 and 0.4. In each case models /lo and I^ 
provide the best agreement with the experimental chord- 
distributions. The numerical results are shown in Ta- 
ble VI and cross-sections of each model shown in Fig. W». 

— (2') (2\ — 

We have plotted Pf-^{r), pexp{r) and measurements of 
the function from the reconstructed samples in Fig. n2. 
The measured data shows some deviation from pr^^{r) 
for p = 0.3. This is due to the accumulation of errors as 
we form the intersection sets of progressively mo re p hase 
functions. Conductivity data is given in Table VII and 
plotted in Fig. |l^. Models I5 and Iio provide a progres- 
sively better estimate of the conductivity. We anticipate 
that increasing the order of model I„ would yield better 
estimates. The results indicate that we have successfully 
reconstructed the 10 S model. 

In Fig. O we have also plotted other data for the lOS 
model. Kim and Torquato |^^ (KT) estimated a for the 
lOS model using a random walker algorithm specifically 
designed to handle locally spherical boundaries. In the 
worst case p = 0.1 our data underestimate that of KT by 
a factor of 1.6 (the error decreases significantly at higher 
volume fractions). This is probably due to the discreti- 
sation effects of our finite difference scheme [g| . 
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FIG. 12. The correlation functions of the lOS model com- 
pared with the "best-fit" function associated with each recon- 
struction. Measurements of p"' obtained from realisations of 
the models are also shown. 

TABLE VII. The conductivity of the lOS model and 
various reconstructions. Models which match the lOS 
chord-distributions (is, 10) provide better estimates of the 
conductivity than a reconstruction based on the single 
level-cut model (Rec. N). 



p 


lOS (KT)=' 


lOS 


Rec. N 


Rec. I5 


Rec. Iio 


0.1 


0.022 


0.014 


0.003 


0.007 


0.011 


0.2 


0.076 


0.063 


0.038 


0.042 


0.052 


0.3 


0.16 


0.14 


0.094 


0.120 


0.13 


0.4 


0.25 


0.24 


0.180 


0.210 


0.22 



^Kim and Torquato, Ref. |^ 

This does not alter our conclusions as all the data pre- 
sented at a given volume fraction are presumably effected 
in the same manner. The data of Bentz and Martys |^ 
(BM) for the lOS model and their one-cut reconstruction 
are consistently lower than ours. 



V. CONCLUSION 





M 


^ 


¥* 

w'-* 


(a) p=0.1 


{b) p=0.2 (c) p=0.3 (d) p=a4 


f- 


m 


m 


m 



(e)p=0.1 {i)p=0.2 (g)p=0.3 (h) ^=0.4 

FIG. 11. The lOS model (a-d) and reconstructions (e-h) 
which reproduce the correlation function (Fig. hd) and 
chord-distributions (Fig. Bft of the model. The conducting 
pore space is shown in black and the images are 10 x 10/xm. 



We have developed a method of reconstructing three- 
dimensional two-phase composite materials from infor- 
mation which can be obtained from digitised micro- 
graphs. First a range of models are generated which 
share low order (volume fraction and two-point correla- 
tion function) statistical properties with the experimen- 
tal sample. The model which most closely reproduces 
the chord-distributions of the experimental material is 
chosen. The distribution functions provided a better sig- 
nature of microstructure than the three-point correlation 
function and are simpler to measure than the microstruc- 
ture parameters C and rj. Significantly the three-point 
and higher order correlation functions of the reconstruc- 
tions can be calculated and employed in rigourous ana- 
lytical microstructure-property relationships. 
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FIG. 13. Conductivity of the lOS model (solid sym- 
bols) compared with various reconstructions (open symbols). 
Model lio provides a very good prediction of the actual con- 
ductivity. Other data are from Refs. @ (KT) and § (BM). 



Three-dimensional realisations of the models can also be 
simply generated for the purpose of numerically evaluat- 
ing macroscopic properties. 

We found that materials with practically identical 
two point correlation functions can have very different 
morphologies and macroscopic properties. This demon- 
strates that reconstructions based on this information 
alone |3P,|^,|l6|-pi| do not necessarily provide a useful 
model of the original material. If the correlation function 
exhibits strong oscillations we found evidence that prior 
methods will provide satisfactory reconstructions. In this 
case it is important to compare the chord-distributions 
of the model and experimental materials. 

Our method can be applied to a wider range of com- 
posite and porous media than prior reconstruction tech- 
niques. The generality of the method is achieved by in- 
corporating new models based on the intersection and 
union sets of level-cut GRF models. The former have 
recently been shown to be applicable to organic aero- 
gels [|3l and porous sandstones |g^ , while the latter may 
be useful for modelling closed-cell foams. Techniques 
based on the single-cut GRF model cannot reproduce 
the low percolation thresholds of these materials p2[. 



The method was successfully used to reconstruct several 
test composites and the overlapping sphere model over a 
range of volume fractions. The reconstructions are better 
able to model the morphology and transport properties 
of the lOS model than prior studies H. 

There are several problems with the reconstruction 
procedure. First, it is possible that two materials with 
different properties may share first and second order sta- 
tistical information and chord-distribution functions. In 
this case the reconstruction method could fail to yield 
good estimates of the macroscopic properties. Second, 
the generality of the models we have employed is not suf- 
ficient to mimic all real composites (although prior stud- 
ies have shown them to be appropriate for a wide range of 
materials P^-pTj). An example is provided above where 
our nine basic reconstructions were unable to model the 
chord-distribution of the lOS model. In this case a fur- 
ther generalisation was found to be successful. Others are 
possible. For example, the restriction that the level-cut 
and length scales parameters are identical for each com- 
ponent of the intersection and union sets can be relaxed, 
or overlapping spheres can be incorporated in the level- 
cut scheme. However the problem remains. It is unlikely, 
for example, that the morphology of randomly packed 
hard spheres could be mimicked by this scheme. Third, 
models formed from the union and intersection sets con- 
tain sharp edges which are energetically unfavourable in 
many materials. However there is little evidence that 
these play a strong role in determining macroscopic prop- 
erties. 

New techniques of characterising microstructure are 
currently being developed such as those based on 
information-entropy pq |. These may contribute to the 
problem of selecting the best reconstructions. Our work 
also has application to the inverse-problem of small-angle 
X-ray scattering from amorphous materials. In this case 
the problem is made more difficult by the absence of 
higher-order information such as chord-distributions (al- 
though some progress may be possible Q ) . Work is un- 
derway to model anisotropic composites and apply the 
method to experimental systems. 
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